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An analog of the Meir-Wingreen formula for the steady-state heat current through a model molec- 
ular junction is derived. The expression relates the heat current to correlation functions of operators 
acting only on the degrees of freedom of the molecular junction. As a result, the macroscopic heat 
reservoirs are not treated explicitly. This allows one to exploit methods based on a reduced descrip- 
tion of the dynamics of a relatively small part of the overall system to evaluate the heat current 
through a molecular junction. The derived expression is applied to calculate the steady-state heat 
current in the weak coupling limit, where Redfield theory is used to describe the reduced dynamics 
of the molecular junction. The results are compared with those of previously developed approximate 
and numerically exact methods. 



T3 

C 

O 

i o i 

(N 
> 

\o 
m 

(N 

o 
o 



X 



I. INTRODUCTION 

Heat transport in nanoscale molecular junctions, i.e., 
in molecules that interconnect metal or semiconductor 
electrodes, is a process that is crucial for the stability of 
the junction and thus for potential molecular electronic 
devices^ — It has been demonstrated experimentally that 
the localized Joule heating may induce a substantial tem- 
perature increase within a molecule-metal contact due to 
inefficient heat dissipation^ Theoretical studies of heat 
transport at the nanoscale, in particular the dependence 
of heat dissipation on various physical parameters of a 
molecular junction, will thus provide valuable insight into 
the transport mechanisms thus facilitating the interpre- 
tation of experimental results and the design of novel 
nanoscale electronic devices. 

Segal and Nitzan investigated the characteristics of 
heat transport of a spin-boson nanojunction model 
(SBNM), where a two- level system is simultaneously 
connected to two heat reservoirs (baths) of different 
temperatures.— ~— Based on these studies, they sug- 
gested possible realization of novel nano-devices such 
as a thermal rectifier^ and a molecular heat pump^ 
The methodology used in their work assumes weak cou- 
pling between the two-level system and the reservoirs so 
that Redfield theor y 11 ' 12 may be applied. This assump- 
tion may not always hold for a realistic molecular junc- 
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tion where the energy flow may be enhanced/maximized 
by tuning certain physical parameters. To address 
this problem, we have developed a numerically exact 
methodology^ to study the dynamics of the SBNM. The 
methodology is based on the multilayer multiconfigu- 
ration time-dependent Hartree (ML-MCTDH) theory,— 
which is a non-perturbative and numerically exact ap- 
proach. Thus, it can provide accurate results in a broader 
range of physical regimes (within the model). For exam- 
ple, our previous study has revealed a turnover behavior 
of the heat current with respect to the coupling strength 
between the two-level system and the heat baths. As a 
consequence, the optimization of heat transport is pos- 
sible by choosing an appropriate set of physical parame- 
ters. 

Numerically exact simulations^ also provide bench- 
mark results that can be used to develop more accurate 
approximate theories. This is the focus of the present 
paper. From a broader perspective, the spin-boson nano- 
junction model is a two-bath, nonequilibrium version of 
the standard spin-boson Hamiltonian for studying elec- 
tron transfer reactions i 15 ' 16 To date, the reduced dynam- 
ics of the two-level system has been of primary interest 
and various approximate approaches have been devel- 
oped to study this dynamics. Examples include the non- 
interacting blip approximation ) 15 ' 17 Redfield theor y 11 ' 12 
Zusman equation^ and its generalization to the low tem- 
perature domain^ It is desirable to apply these methods 
directly to the study of heat transport in the spin-boson 
nanojunction model. However, this is not straightfor- 
ward since the operator for heat current involves not only 
the degrees of freedom of the two-level system but also 
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those of the heat baths. 

On the other hand, the relation between the reduced 
dynamics of an open system and the current flowing 
through it was found in the theory of electron transport 
in mesoscopic systems by Meir and Wingreen.— They 
derived a Landauer-type expression that relates the elec- 
trical current to certain correlation functions of the small 
mesoscopic system. The goal of this paper is to derive 
a similar expression for the heat current through a two- 
level system, driven out of equilibrium by two heat baths. 
Based on the thus obtained expression, the approximate 
methods mentioned above can be used to study heat 
transport properties of the SBNM. 

The remainder of the paper is organized as follows. 
The model Hamiltonian and the observable of interest 
are described in Sec. Ill Al Section III Bl describes the the- 
oretical techniques that allow us to express the heat cur- 
rent through correlation functions of the two-level sys- 
tem. The final Meir- Wingreen type expression for the 
steady-state heat current is obtained in Section III CI 
As a demonstration, we apply the developed theory in 
Scc. lIIDl to evaluate the heat current in the limit of weak 
coupling between the two-level system and the baths. In 
Sec. IIIII this heat current is compared with the results 
obtained using the Segal-Nitzan approach and the nu- 
merically exact simulations. Sec IIVI concludes. 



II. THEORY 



A. Spin-Boson nanojunction model 



paper) 



Hb — He + Hjj , 
H K = ^2 u m al H a m ; K = C,H, (2) 



me A' 



where (a m ) are bosonic creation (annihilation) oper- 
ators. The bridge Hamiltonian H$ in the second quanti- 
zation reads 



i=l,2 



where riy = c^o, and cj (cj) is a fermionic creation (anni- 
hilation) operator. The energy spacing between the two 
levels of the bridge is denoted by e. The coupling between 
the bath and the TLS is given by 

Hsb = Hsc + Hsh, 
H S k = V m(<4n + a m )[n 2 i + n 12 \; K = C, H. (4) 



EiTiu) E 2 - E 1 = e > 0, 



(3) 



The dynamics of the TLS is restricted to the single- 
particle space spanned by the basis functions |1) = 
cjjvac) and |2) = <4|vac), where |vac) denotes the 
fermionic vacuum. This restriction guarantees that the 
TLS representation by Hamiltonian ([3]) is equivalent to 
the more commonly used form ) 15 ' 16 Hs = |(|2)(2| — 
1}(1|). On the other hand, the fermionic representation 
introduced here is more convenient for the diagrammatic 
technique exploited in Sec. Ill Bl 



The spin-boson nanojunction model (SBNM) consid- 
ered in this worke r 10 ' 13 consists of two heat reservoirs 
(baths) interconnected by a bridge, which is represented 
by a two-level system (TLS) (see Fig.Q]). The heat baths 



Starting with the Heisenberg operator of the heat cur- 
rent from the TLS to the cold bath 



Ic(t) = i[H,H c }(t) 



(*), (5) 
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FIG. 1: Spin-Boson nanojunction model. 



at different temperatures drive the TLS out of equilib- 
rium and the dependence of the heat current through the 
TLS on various parameters of the model can be studied. 
The SBNM Hamiltonian reads as 



H = Hb + Hs + Hsb , 



(1) 



where Hb describes the two harmonic baths, "cold" (C) 
and "hot" (H) (atomic units are assumed throughout the 



straightforward algebra gives the expectation value of the 
heat current 



(l c (t)) = -21m 



i.j; m£C 



(6) 



In this expression we introduce the notation 

C = (i-<yv m , (?) 

and the nonequilibrium Green's function (NEGF) 
M mt ij(t,t') defined on the Schwinger-Keldysh 
contour— iH (depicted in Fig. [5]) 



M m ,ij (*,*') = {T c [a m (t)n tJ {t')\) 



(8a) 



M^-ftf); t> c t>, 
M< y (t,f); t<J 
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where t > c t' means that t is "later" on the contour 
than t' as is illustrated in Fig. [5] The contour-ordering 
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FIG. 2: The Schwinger-Keldysh contour. 



operator T c moves operators with later (in the contour 
sense) time-arguments to the left. Thus, in accordance 
with Eq. ©, the greater and the lesser Green's func- 
tions are dchned as M^ i -(t,i / ) = {a m (t)riij(t')) and 
^mij = ( n ij{t') a m(t)), respectively. The choice 
of t max , the contour's turning point, is arbitrary as long 
as t max > t,t'. The quantum mechanical average with 
respect to the initial state in Eq. (|8a|) is defined as 
(...) = Tr[... p(to)}, where the initial state density op- 
erator p(to) is taken in the product form 



p(to) 



PB X p S , 



PB 



e ~PcH c e -0 H H H 



3 c H Cp -PhH 



Tr {e- 



(9a) 



(9b) 



In the expression above ftc = V^B^c, Ph = ^-/ksTn, 
and ks is the Boltzmann constant. The initial density 
operator of the bridge ps is arbitrary as long as one is 
only interested in the steady state. 



B. Dyson equation for AI m ,ij(t,t') 



Inspired by the pioneering work of Mcir and 
Wingreen,— we proceed to express M mi ij(t,t') in terms 
of correlation functions of the TLS operators only. Our 
derivation closely follows the procedure used by Gaudin 
in his proof of the generalized Wick's theorem ! 23 ' 24 First, 
we rewrite M m ^j(t,t') in the interaction representation 



M mtij {t,t') = Tr{S c [a m (t)n i:j (t')]p(to)} ■ 



(10) 



Operators in the interaction representation are denoted 
by a hat, i.e., o m (r) = e i{ - Hs+H ^ 7 'a m e~^ Hs+H ^ T , and 
S c is 



1 + H) / dTH SB {T) 



dr / dr' H S b(t)Hsb(t) 



(11) 



The above expansion of S c in powers of Hsb(t) yields 
the perturbation expansion for M m ^j{t,t') 

M mtij (t, t') = M£% (t, t') + Mfy (t, t') + ■ ■ • , (12) 
where 

M^t') = Tr^a^n^t'Mto)} = 0, 
M^(t,t') = -i J dr Tr[T c [H S B{T)d m (t)n l0 {t')]p(t Q )} , 

(13) 

The successive commutation of a m through all operators 
Hsb and p(to), with the use of the following identit y 23 ' 24 

exp(-/3o; m aJ rl a m a m ) = a m cxp(-/3w m a] n a m ) exp(-/3w m ), 

(14) 

and the cyclic permutation within the trace, leads to 



M^(t,t') = -i^V^ / dr Tr{T c [a m (t)al(r)]p(t )} 
ki Jc 
xTr {T c [n k i(T)n l3 (t')]p(t )} , 

M%Ut,t>) = (-i) 2 J2 V ™ I drTr{T c [a m (t)al(r)]p(to)} 



kl 



x / dr 



' Tr{T c [# SB (rO^)M*')]/>(*o)} , 



(15) 



By continuing the series in Eq. (| 15[) one notices that 

^Sy(M')=E^ [drD m (t,r)K%P(r,t'), (16) 
ki Jc 

where 

D m {t,r) = -iTr {T c {d m (t)al(T)}p(t a )} (17) 

is the contour-ordered Green's function of the non- 
interacting bath and ;y(T, t') arises from the pertur- 
bation expansion of the Green's function 

oo 

K M ,ij(r,f) = J2 K $ij( T >^ = ^{T c [nki{T)n l3 (t')] P {tv)} 

q=0 

(18) 

Finally, substitution of Eqs. (02) and ([HJ into Eq. ([12]) 
results in 

M mtij (t,t') =J2 V m f dr D m (t,T)K kltij (T,t'), (19) 
ki Jc 

which expresses M m ij(t,t') in terms of correlation func- 
tions of TLS operators and non-interacting bath oper- 
ators separately. It is noted that Eq. (jT9")) can also be 
obtained using the equation-of-motion technique.— ~— 

A Dyson equation similar to Eq. (|19[) was recently ob- 
tained in the context of charge transport through sin- 
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gle molecules, where molecular many-body states were 
employed to describe the central molecular part^ This 
similarity is due to the formal resemblance of the Hamil- 
tonian describing the bath-TLS coupling, Eq. ^ in the 
present work, and that of molecule-contact coupling, Eq. 
(8) in Ref. S 



transform we finally obtain 



id 



(^) = E E "rnV%V m 
i,j;k,l m£C 



0- + r&)Ku,ijM 



(24) 



where 



C. Steady-state heat current 

Since only the greater NEGF M>^-(f, 4') is needed 
to evaluate the heat current in Eq. ©, we can im- 
mediately assign t and t' to the lower and the up- 
per Schwinger-Kcldysh contour branches, respectively 
Omitting straightforward manipulations with Eq. (|19[) 
based on Langreth's contour deformation rules^ the fi- 
nal result for Al> ..(f, f'). which contains only traditional 
real-time correlation functions and the time integration 
along the real axis, is obtained 



M>iAt,t')= > -v. 



^ykl 

kl 



dr 



la 



D r m (t,T)K>(r,t>) 



+D>(t,T)K% l>ij (T,t>)] 



(20) 



Here, correlation functions D(t,r) and K(T,t') are given 
in Eqs. (|17[) and (|18[) . respectively, and for an arbitrary 
contour-ordered Green's function G(t,t') we define 

G r (t,t') = 6(t-t')[G > (t,t')-G < (t,t')}, 

G a (t, t') = 8{t' - t) [G< (t, t')-G> (t, t')] , (21) 

where 8(t — t') = 1 if t > t' and zero otherwise. 

In the steady-state regime (to — > — oo) two-time cor- 
relation functions become stationary, i.e., M>^(^,i') = 
A/„; ,.;/ - t 1 ). Combining Eq. @ with Eq. ^ and set- 
ting t = for simplicity wc obtain the expression for the 
heat current 



(Ic) 



- 2Im E 



u m v*>v*? 



dr (D' r m (-T)K> t3 {T) + Dl{-r)KZ Uij {T)) . 

(22) 



The Green's functions of the cold non-interacting bath 
can be easily evaluated a o 26 ' 29 



D r m (-r) = -^(-r)e-^, 
D>(-r) = -i(l + r ] g)e i ^, 



(23) 



where 77^ = ri c (uj m ) ~ [e /3cLJm — l]" 1 is the Bose-Einstein 
distribution. Combining the last two equations and not- 
ing that the time integration corresponds to a Fourier 



+00 



(tK 



(25) 



The expression for the heat current to the other (hot) 
bath is obtained by a straightforward substitution of "C" 
with "H" in Eq. (JMJ) 

Eq. (|24p exactly relates the heat current in the spin- 
boson nanojunction model to correlation functions of the 
bridge that can be evaluated using methods developed 
to describe reduced dynamics of a TLS. In fact, Eq. ([2~4")l 
is also formally valid for a multi-level bridge. In this 
case the indices i,j,k,l run through all multiple bridge 
states and the expression includes the corresponding cou- 
pling constants Vjfi and correlation functions K^^(lu). 
A similar expression was recently derived for photonic 
heat current through an arbitrary (nonlinear) circuit 
element coupled to two dissipative reservoirs at finite 
temperaturesi 30 i 31 

The next subsection is dedicated to an approximate 
evaluation of the TLS correlation functions, and hence 
the SBNM heat current, in the limit of weak TLS-bath 
coupling. 



D. Treatment of the correlation function within 
the Redfield approximation 

Redfield theory in the form of kinetic equations for the 
TLS populations has been previously employed by Se- 
gal and Nitzan to study the heat-conducting properties 
of the SBNM^i^ Within the Redfield approximation it 
is possible to obtain simple and physically transparent 
analytical results. However, this approach is accurate 
only at very weak TLS-bath couplings, which restricts 
its applicability. In particular, the non-monotonic de- 
pendence of the heat current on the coupling strength 
(see the next section for details) can not be described 
within this approach^ 

In the remainder of the paper we will use the theory 
developed in the previous subsections to derive an im- 
proved description of heat transport that still retains the 
analytical nature of the method. Specifically, we will 
evaluate the TLS correlation functions K^^(lu) within 
the Redfield approximation. The thus evaluated corre- 
lation functions will then be used to calculate the heat 
current by means of Eq. (f24|) . As will be demonstrated 
in Sec. Mil this approach results in a notable improve- 
ment over the original Rcdficld-typc theory by Segal and 
Nitzan. 
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In the following discussions it is convenient to split with 



(t,f) into the "retarded" (+) and "advanced" (— ) 

part 

K> Aj (t,t') = K+.^t') + K-^t'), (26) 



where 
K 



(27a) 

Aj (t,t') = e(t'-t)Tr{e iH ^-^n ij e- iH ( t '-^p(t)n kl } . 

(27b) 

Lesser correlation function K kl i At, t') can always be ob- 
tained using the property K kl . .(t, t') 



K 



K^'ty Fur - 



thermore, the retarded and the advanced parts are re- 



K + 

H,ij 



(t,t>) 



so we focus 



lated through K ji[k (t',t) 
on the retarded correlation function and obtain all other 
functions accordingly. The steady-state regime implies 
that both t — to and t' — to significantly exceed the char- 
acteristic time of the transient dynamics. This, together 
with the Redficld approximation ) 11 ' 12 yields the steady- 
state density operators in Eq. (|27|) 

p(t) = p(t<) = (|l>(l|i? + \2}{2\P 2 S ) p B , (28) 

where P[ , P§ are the stationary populations of the two 
levels and ps is defined in Eq (|9bl) . Since the density 
operators p(t) and p(t') in Eq. (|2"5|) arc effectively time- 
independent, the correlation functions in Eq. (|27|) depend 
only on r = t — t' as expected at the steady state. 



The correlation functions kl have a clear physi- 
cal meaning. For example K^ 2 12 (t) describes a process 
where the system, being in the stationary state, under- 
goes a coherence- introducing transition |2)(2| — > |1)(2|, 
then evolves for time t, and finally the amplitude of 
the coherence |2}(1| is measured and multiplied by the 
steady-state population of state |2}(2|, i.e., P 2 ■ There- 
fore, within the Redficld approximation, this correlation 
function can be calculated as 



K+(T) = Pi'(2\p s (T)\l), 



(29) 



where ps(i~) is the solution of the Redficld master equa- 
tion 



dr 



= -iEijPij(T) + 2J RijkiPki{r) (30) 



hi 



subject to the initial condition ps(0) = |1)(2|. In this 
master equation, pij — {i\p{T)\j}, Eij — Ei~ Ej and the 
Redfield tensor Rijki is given by 

Rijkl = ^tjih + T ijih - S M r ^"fc " Slk T 7aaj' ( 31 ) 



r tjik = / d S Tr B {(l\e- iH ° s H SB e iH ° s \j)(i\Hs B \k)} 
Jo 



-iE ik s 



/>oo 

r,« fc = / dsTi B {(l\H S B\j)(i\e- iHsB Hs B e iHBS \k)} 
Jo 

(32) 



x e 



-iEijS 



All the other retarded correlation functions can be ob- 
tained in a similar fashion. 



Instead of solving Eq. (|30|) for each correlation func- 
tion independently, it is more practical to reformulate 
the problem as an equation of motion for the correlation 
functions themselves. Straightforward but somewhat te- 
dious algebraic manipulations, including the evaluation 
of various components T^ fc -* of the Redfield tensor, al- 
low one to obtain this equation of motion in matrix form 
with the initial conditions included 

#K+(r)=f ie '~2 iAe , + 7 V(r) 
dr K-iAe + j -ie' - 7 / v ' 



(1 Pi 



s > 



where K + is 2 x 2 matrix 



(33) 



(34) 



The renormalized energy splitting e' of the TLS levels is 
given by 

e' = e + Ae, 

1 f°° 
Ae=- V V duj[l + 2r] K {uj)} 

77 K tt,H J o 



X Jif(w) 



e + ui e — u> 



(35) 



where V stands for the Cauchy principal value. The spec- 
tral density 

J K (u) = 7T V£5(u - w m ); K = C,H, (36) 

is introduced here to change the summation over the bath 
modes to an integration over the frequency in Eq. (|35l) . 
In this work we employ a spectral density of Ohmic form 
with exponential cutoff 



Jk(u) = \k~. e" 



(37) 



where and u>k are the TLS-bath coupling strength 
and the characteristic frequency of the K th bath, respec- 
tively The stationary dephasing rate is related to the 
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spectral density as 



7= [1 + 2^(6)1^(6) 

K=C,H 



(38) 



The system of linear ordinary differential equations in 
Eq. (|33|) can be easily solved by converting it to a system 
of linear algebraic equations employing a Fourier trans- 
form 



f 



D(u) 



~(~f + ie'-iuj)ps, ( 7 + iAe)Pf 
(7 - iAe)Pg, (7 - ie' - iu)Pf 

(39) 

where D(u) = (e') 2 — (Ae) 2 — 2iuj — uj 2 . The inverse 
Fourier transform is not needed since these correlation 
functions in frequency domain are exactly what is re- 
quired to evaluate the heat current in Eq. (|24|) . The 
expression for K + is the only one explicitly needed since 
K = (K + )^ and the greater and lesser correlation func- 
tions are given by 



K 



kl.ij 



(40) 



Finally, Eq. (|2"4")l can be rewritten with the use of 
Eqs. @g) and p}]) 



(^C) = - d . 2 

where heat transfer rates are 



(41) 



r c - 



r° = -- 

7T 



Jc(w)wxRe 



dwJc(w)wxRe 



27 — ie(l + 2rjc(ui)) — iu> 
(42a) 

2 7 -ne(l + 2T7 C (u)) -iu 
D(w) 



(42b) 

The heat transfer rates for the hot bath can be found 
from the expressions above by formal substitution of "C" 
with "H" . 

In a weak coupling regime the heat transfer rates re- 
duce to 



T c d =2eJ c (e)[l + V c(e)}, T, c 



2eJ c (e)r,c(t), (43) 



which arc identical to the results of Ref. y. Another way 
to obtain the rates in Eq. (|43|) is to eliminate the de- 
phasing (and energy rcnormalization) by setting 7 and 
Ae in Eq. (|33[) to zero. Hence, the previous expression^ 
for the heat transfer rate is recovered when coherence, 
described by correlation functions K(t), persists indef- 
initely. Once realistic dephasing and energy renormal- 
ization are added, Eq. (|4"2")l for the heat current is ob- 
tained. From this we expect that Eq. (|42p can provide 
a more accurate description of the SBNM heat current 
than Eq. (|43p except for the very weak coupling regime 
where they agree with each other. In what follows, we 
will refer to the expression for the heat current obtained 



in this work [Eqs. gj), ([4"2])] as NEGF-Redficld approach 
to distinguish it from the Rcdfield approximation used 
previously^ 

Finally, it is worthwhile to discuss the energy conserva- 
tion in the NEGF-Redficld approach. Conservation laws 
are not always automatically satisfied in approximate 
theories. Thus, often special care must be taken to guar- 
antee conservation of, e.g., number of particles, energy 
or momentum. For example, in the application of many- 
body Green's functions it has been demonstrated that 
self-consistency must be incorporated into the Dyson 
equation in order to preserve conservation laws . 32 ' 33 As 
will be seen later, the NEGF-Redficld approach devel- 
oped in this subsection does not conserve energy, i.e., 
(Ic) + (Ih) 7^ 0, if the TLS steady-state populations are 
determined as the stationary solution of Eq. (|3"U)) . In fact, 
numerical tests show that with so obtained steady-state 
populations the heat current in Eq. (|4"Tj) does not nec- 
essarily vanish at equilibrium, i.e., when the both baths 
have identical temperature. To correct for this, we intro- 
duce self-consistency by explicitly requiring the energy 
conservation 



(Ic) + (Ih) = 0, 



(44) 



The steady-state populations are then determined from 
this constraint. Specifically, once the heat transfer rates, 
Eq. (|4"2")l . arc evaluated for both baths, Eq. (|4"4")l becomes 
a linear equation for the steady-state populations of the 
two levels. This linear equation is solved straightfor- 
wardly taking into account the identity Pf + P§ = 1. It 
is noted that the so obtained steady-state populations are 
identical to those obtained by assuming a kinetic equa- 
tion for the TLS populations with the "cold" and "hot" 
rates given by Eq. (|4"2")l . The latter procedure can be 
straightforwardly extended to multi-level bridges. 

It will be demonstrated in the next section that the in- 
troduction of self-consistency performs surprisingly well 
and drastically improves the accuracy of the NEGF- 
Rcdfield method. 



III. NUMERICAL RESULTS AND DISCUSSION 

The dependence of the steady-state heat current (Ic) 
on the TLS-bath coupling strength Ac is shown in Fig- 
ure G3 The bath parameters in Eq. (|37|) . i.e., the TLS- 
bath coupling strength and the characteristic frequency, 
are chosen to be identical for the two baths, i.e., Ac = Xh 
and ix>c = Wjj, respectively. The bath temperatures are 
defined by k B T c /u} C = 0.1, k B T H /u} H = 0.21 for the 
cold and the hot bath in panel (a) and fc^Tc/wc = 0.17, 
ksTH/^H = 0.26 for the cold and the hot bath in panel 
(b), respectively. The ratio of the TLS energy spacing 
to the characteristic frequency of the bath is taken as 
e/ujc = 0.75 in both panels. 

The numerical results of the NEGF-Redficld and Segal- 
Nitzan methods are obtained from Eq. (l4"2"T) and Eq. (|43p , 
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FIG. 3: Dependence of the steady-state heat current Ic on 
TLS-bath coupling strength Ac- Panels (a) and (b) corre- 
spond to different temperatures of baths. The spectral den- 
sities for the two baths are equal (Ac = Ah, wc = wsr)- The 
TLS energy spacing is e/tuc = 0.75. The numerical results 
have been obtained with the numerically exact ML-MCTDH 
methodology (circles), the Segal-Nitzan theory (dashed lines), 
and the NEGF-Redfield approach developed in this work (full 
lines). Error bars for some numerically exact points are 
smaller than the size of circles, and, therefore, not shown in 
the figure. 



respectively. The numerically exact results are obtained 
by means of the Multilayer Multiconfiguration Time- 
Dependent Hartree (ML-MCTDH) approach^ The de- 
tails on the application of the ML-MCTDH method to 
the study of heat transport in the SBNM are described 
elsewhere For the parameters considered here both the 
NEGF-Rcdficld and Segal-Nitzan theories, depicted by 
the dashed and full black lines, perform well if Ac is suffi- 
ciently small. However, when Ac/wc is larger than ~ 0.5, 
both approximate theories cease to be valid and the re- 
sults deviate significantly from the numerically exact sim- 
ulation. This failure is expected considering the weak- 
coupling nature of the approximations in both approxi- 
mate methods. In the intermediate region (Xc /loc ~ 0.3- 
0.5) the NEGF-Redfield method gives significantly better 
results than the Segal-Nitzan approach, because the for- 
mer treats the TLS-bath coupling more accurately as was 
discussed in the previous section. 

At large coupling strengths the numerically exact sim- 
ulation predicts a "turnover" t 13 i 34 i 35 i.e., the heat cur- 
rent reaches its maximum and starts to decrease with 
increasing of the coupling strength. This phenomenon 
is similar to Kramers' turnover and has been discussed 
in more detail elsewhere^ The turnover is not repro- 
duced (even qualitatively) by the approximate theories 
although the NEGF-Redfield method displays the correct 
curvature change in the intermediate coupling regime. 

The accuracy of the NEGF-Rcdficld approach dete- 
riorates significantly if the self-consistency constraint, 



Eq. (|44p. is not applied and the steady-state popula- 
tions are obtained from regular Redfield theory, i.e., from 
Eq. (|3"D"|) . The heat currents from the bridge to the 
cold bath, (Ic), and from the hot bath to the bridge, 
— (Ih), evaluated for this choice of the steady-state popu- 
lations, are depicted by up-triangles and down-trianglcs, 
respectively, in Fig. H^a). In the inset, it is seen that 
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FIG. 4: Dependence of the steady-state heat current (Panel a) 
and the steady-state population of the higher TLS level (Panel 
b) on the TLS-bath coupling strength Ac- All parameters are 
the same and circles, full lines and dashed lines correspond to 
the same theoretical methods as in Fig. [2 a). The up-triangles 
and down-triangles depict (Ic) and —{Ih}, respectively, in 
the case when steady-state populations are obtained directly 
from the Redfield master equation, Eq. (|30|l , and not from the 
self-consistency requirement, Eq. (|44[) , The average of these 
two currents, i.e., ({Ic} — (7h))/2, is depicted by the dashed- 
dotted line. The thin red line in Panel b, which was obtained 
by fitting the numerically exact results, is a guide for the eye. 



at very small coupling strength these two currents co- 
incide with each other and with the current determined 
self-consistently. However, {Ic) and —(Ih) start to de- 
viate from each other at Xc/toc ~ 0.01, thus demon- 
strating the non-conservation of the energy current dis- 
cussed above. Furthermore, the heat current from the 
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hot bath to the bridge becomes vmphysically negative at 
\ c /uj c ~ 0.1. 

The energy conservation condition in Eq. (pHf suggests 
another, more symmetric, definition of heat current by 
averaging (I c ) and ~{I H ), i.e., {(I c ) - (I H ))/2M This 
averaged current, with (Ic) and —(Ih) given by up- and 
down-triangles, respectively, is depicted by the dashed- 
dotted line in Fig. 2J This current is seen to agree much 
better with the numerically exact results than (Ic) and 
— (Ih) separately. However, this "averaging" trick does 
not resolve the energy non-conservation problem, and 
merely conceals it. Furthermore, the NEGF-Redfield ap- 
proach with self-consistency introduced (full line) is still 
in better agreement with the numerically exact calcula- 
tions. Therefore, the results of Fig. HJa) emphasize the 
necessity of the energy conservation requirement and sup- 
port our choice of the steady-state populations. 

This choice is further supported by Fig.QJb), where the 
steady-state population of the higher TLS state, eval- 
uated from Eq. flM]) within the NEGF-Redfield theory 
(full line), Redfield master equation (dashed line) and 
the numerically exact simulation (circles), is depicted. 
Whereas standard Redfield theory (dashed line) yields a 
population that is totally independent on the cou- 
pling strength, the NEGF-Redfield theory performs bet- 
ter reproducing, although not quantitatively, the increase 
of the population with Ac- 

The influence of the TLS energy spacing e on the 
steady-state heat current in shown in Figure [S] As is 
seen, the heat current also exhibits a turnover behavior 
as a function of the energy spacing. This is due to the 
resonant character of heat transport: The heat transport 
is most efficient if the TLS has the same characteristic 
frequency, i.e., energy spacing between the two levels, as 
that of the baths. Too large or too small e (compared to 
w c ) both result in inefficient heat transfer. The detailed 
discussion of this phenomenon is given elsewhere 

While the Segal-Nitzan theory can qualitatively pre- 
dict this turnover behavior, there is a significant quan- 
titative discrepancy with the numerically exact simula- 
tion. However, the NEGF-Redfield theory method is in 
much better agreement with the numerically exact results 
in Fig. [5] We attribute this to the fact that since the 
coupling strength is relatively small (Ac/wc = 0.25), a 
more accurate treatment of TLS-bath coupling within the 
NEGF-Redfield approach leads to a much better agree- 
ment with the numerically exact result. 



IV. CONCLUDING REMARKS 

In this paper we have derived an analog of the Mcir- 
Wingreen formula for a spin-boson nanojunction model 
for heat transport in a single-molecule junction. This 
formula exactly relates the steady-state heat current 
through a two-level system, connected to two heat reser- 
voirs, to correlation functions of the operators of the two- 
level system only. The formula allows one to calculate the 
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FIG. 5: Dependence of steady-state heat current Ic on the 
TLS energy spacing e. The results have been obtained with 
the numerically exact ML-MCTDH method (circles), the 
Segal-Nitzan theory (dashed line), and the NEGF-Redfield 
approach (full line). The temperatures are feTb/wc = 0.17 
and kBTn/ojH = 0.26 (ljc — ^h)- The TLS-bath coupling 
strength is Xc/uc = 0.25 (Ac = A_h). Error bars for some 
numerically exact points are smaller than the size of circles, 
and, therefore, not shown in the figure. 



steady-state heat current for the spin-boson nanojunc- 
tion model using previously developed methodologies for 
the reduced dynamics of the standard spin-boson prob- 
lem. 

As an illustrative application of the formalism devel- 
oped, we have analyzed the properties of the steady-state 
heat conduction in the spin-boson nanojunction model 
using the derived expression for the heat current and 
Redfield theory to evaluate the correlation functions of 
the two-level system. In addition, a self-consistency cri- 
terion was formulated and applied that enforces the con- 
servation of the heat current within the NEGF-Redfield 
scheme. Employing this approach, we have studied the 
dependence of the heat current on the energy spacing 
of the two level-system and the coupling strength be- 
tween the system and the heat baths. The numeri- 
cal results obtained in this work demonstrated that the 
NEGF-Redfield method represents a significant improve- 
ment compared to previous approaches due to the more 
accurate treatment of the coupling between the two-level 
system and the heat baths. 

Finally, we would like to comment on the relation of 
the Meir-Wingreen formula ([2^]) developed in this work 
to previously developed NEGF-based theories of thermal 
transport. The latter theories are based on the phonon 
picture and the anharmonicity, i.e., phonon-phonon in- 
teractions, is included perturbatively (see Ref. [36| for 
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an excellent review and references therein). In con- 
trast, Eq. (f2"4")l treats in-bridge anharmonicity non- 
perturbatively since exact eigenstates of an anharmonic 
oscillator can serve as bridge levels. For example, the 
two levels considered in this may model the two lowest 
eigenstates of a double- well potential, which often cannot 
be treated perturbativcly. On the other hand, an eigen- 
state expansion is limited to a relatively small bridge in- 
cluding only few vibrational modes. Thus, the approach 
developed in this paper has to be extended and opti- 
mized to apply it to large molecular bridges. Therefore, 
the method developed in this work and the conventional 
NEGF-based theories of heat transport are independent 
and complementary. A major improvement in the the- 
ory of heat transport through nano junctions can thus 
be expected if the two methods arc combined, allowing 
the simultaneous treatment of weakly and highly anhar- 
monic modes of the same nano junction by anharmonic 
perturbation theory and exact diagonalization methods, 



respectively. 
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